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The structure and dynamics of the water/vapour interface is revisited by means of path-integral 
and second-generation Car-Parrinello ab-initio molecular dynamics simulations in conjunction with 
an instantaneous surface definition [A. P. Willard and D. Chandler, J. Phys. Chem. B 114, 

1954 (2010)]. In agreement with previous studies, we find that one of the OH bonds of the water 
molecules in the topmost layer is pointing out of the water into the vapor phase, while the orientation 
of the underlying layer is reversed. Therebetween, an additional water layer is detected, where the 
molecules are aligned parallel to the instantaneous water surface. 


INTRODUCTION 


The water/vapor interface is ubiquitous in many in¬ 
teresting applications involving both chemical processes 
and phenomena in biology and aqueous chemistry [i|. A 
better understanding of the interfacial behavior of water 
is crucial for studying phenomena as diverse as protein 
folding [3, Q, the structure and function of biological 
membranes and membrane proteins [^, the Hofmeister 
effect 0, electrochemical processes in aqueous batteries 
0, and the remarkable organic catalysis on water 
to mention just a few. Despite long study, understanding 
the molecular aspects of hydrophobic solvation continues 
to be an area of active research [lol - [l^ . 


The simplest case of an aqueous-hydrophobic bound¬ 
ary is the water/vapor interface, where the vacuum can 
be considered as the “ultimate hydrophobe” [l^ , with 
no attractive or repulsive van der Waals forces between 
the water and the hydrophobic phase. Due to its pro¬ 
totypical character, this particular interface has been 
studied extensively, both experimentally [SSl and 
theoretically 3ll-l54|. Moreover, there is an ongoing de¬ 
bate in the literature regarding the pH of water at the 
water/vapor interface. Different experimental and the¬ 
oretical studies have suggested the possibility that the 
pH at the surface of water could be either basic or acidic 
More recently, it has also been proposed that 


enhanced charge transfer at the interface could rational¬ 


ize the existence of the negative charge at the surface of 
water W, 6 ^. In any case, even the bare water/vapor 
interface remains highly contentious. 


The broken symmetry at the water/vapor interface re¬ 
stricts the ability of water molecules to form H-bonds 
like they do in the bulk. One of the most common 
quantities used to characterize the average configuration 
of water molecules is the ratio of water molecules that 
corresponds to single-donor (SD), double-donor (DD), 
or non-donor (ND) configurations. Within this frame¬ 
work, perfect crystalline bulk ice consists of 100% DD 
configurations, whereas bulk liquid water and its sur¬ 
face are expected to have a mixture of all three con¬ 
figurations with varying ratios induced by thermal fluc¬ 
tuations. Various experimental techniques have been 
applied in order to probe the structural properties at 
the interface. In particular, using surface-sensitive vi¬ 
brational sum-frequency generation (SFG) spectroscopy, 
Shen and coworkers observed a complete suppression of 
the free OH peak at 3700 cm“^ by titrating the dangling 
OH groups at the water/vapor interface with methanol 
ll^ . From these measurements, they estimated that 
the vapor/water interface consists of 25% SD and 75% 
DD configurations. This implies that a significant frac¬ 
tion of the water molecules at the interface possess dan¬ 
gling OH bonds that are protruding out of the water 
phase into the vacuum. In later studies, dangling OH 
groups at the interface have been confirmed by SFG mea- 
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surements of others 22, and molecular dy¬ 

namics (MD) simulations 31. 33. 34. 36l 38. 45. HsL 6^. 
Nevertheless, other experiments using X-ray absorption 
[nm and SFG spectroscopy [i^, as well as ab initio 
MD (AIMD) simulations [1^ have found a much larger 
fraction of “accepter-only” water configurations. How¬ 
ever, quantifying the relative occurrence of the various 
water configurations entails a multitude of computational 
challenges, which are due to the large dipole moment, 
polarizability and nuclear quantum effects (NQE) of liq¬ 
uid water, as well as nonadditive cooperative effects of 
H-bonds Hence, the predictive power of such sim¬ 


ulations critically depends on the accuracy of the em¬ 
ployed interaction potential, but also on the involved 
time and length scales to minimize statistical uncertain¬ 
ties and single-particle finite-size effects [67|. In case of 
the water/vapor interface, an even larger system has to 
be considered in order to stabilize a two-dimensional pe¬ 
riodic water slab with a sufficiently large vacuum portion 
on top of it to eliminate spurious long-range interactions 
at the surface. 

In an earlier work, we studied the water/vapor in¬ 
terface from first principles usi^ the second-generation 
Car-Parrinello AIMD method ^ , 68, |6^ . In this study, 
no evidence of a significant occurrence of ND water con¬ 
figurations at the water surface was observed. However, 
determining the exact location of the water/vapor inter¬ 
face is quite ambiguous due to the inherently large spatial 
and temporal fluctuations of the instantaneous interface. 
Yet, the particular definition of the water surface can sig¬ 
nificantly affect the quantitative values of DD, SD and 
ND water molecules near the water/vapor interface. As 
a consequence, there have been several attempts to refine 
the definition of the interface at a molecular level that ac¬ 
counts for the instantaneous fluctuations [zB-Izl- More 
recently, an instantaneous interface definition to eluci¬ 
date the presence of enhanced structural correlations at 
the interface that are not captured by a criterion based 
on a time-averaged density profile has been proposed [74 1. 

However, since water mainly consists in a large part 
of light H atoms, NQE, such as the quantum mechanical 
zero-point energy (ZPE) and tunneling effects, are po¬ 
tentially essential to obtain the correct quantitative, and 
sometimes even qualitative behavior, and thus should be 
explicitly taken into account l75l482l| . Although it is 
well known that NQE generally weaken intermolecular 
H-bonding, resulting in a faster rotational and transla¬ 
tional dynamics and at the same time less structured 
liquid, there is an on goin g debate regarding the magni¬ 


tude of this effect 82-861. In fact, Habershon et al. have 


shown that there are competing NQE in bulk water: on 
the one hand the enhanced quantum fluctuations of the 
protons result in stronger H-bonding, increased structure 
and slower dynamics, while on the other hand, the en¬ 
hanced librational modes weakens the H-bonds, reduces 
the structure and increases the diffusion |86|. Moreover, 


there exists considerable interest on the impact of NQE 
on the water/vapor interface. In particular, both ex¬ 
perimental and theoretical studies have shed light on the 
presence of isotope fractionation effects observed between 
the liquid water and vapor phase, which is due to NQE 
[87 h 91|| . In general, the role of NQE on the structural 
and dynamical properties of the water/vapor interface 
remains poorly understood. 

In this work, we revisit the orientation of water 
molecules at the water/vapor interface, as well as the 
structure and dynamics of the corresponding H-bond net¬ 
work based on AIMD and including NQE, by means of 
path-integral MD (PIMD) simulations. At variance to 
our previous study [i^ , the improved instantaneous sur¬ 
face definition of Willard and Chandler is employed to 
identify the individual water layers and to investigate 
the orientational distribution, H-bond network, as well 
as dynamics of interfacial water. 

The remainder of this paper is organized as follows. 
In Section 2, we summarize the finite temperature path- 
integral technique to account for ZPE and tunnelling ef¬ 
fects, as well as the H-bond and instantaneous surface 
definitions. Thereafter, in Section 3, we describe the 
computational details of our AIMD and PIMD simula¬ 
tions. Our results are presented and discussed in Section 
4 before concluding the paper in Section 5. 


COMPUTATIONAL METHODS 


Path-Integral Formalism 


In the PIMD method, all quantum particle are re¬ 
placed by harmonic P-bead ring-polymers, which are 
then treated classically. This is to say that quantum me¬ 
chanical properties can be exactly calculated by sampling 
the path-integral phase space using MD, since the ex¬ 
tended ring-polymer system is isomorphic to the original 
quantum system [92l494l |. For this purpose, the canonical 
quantum partition function Z[j3) is expressed in terms of 
the inverse temperature = /c_bF, i.e. 


(e-/3pA)^ 


= lim Zp(^), 
P—fOO 


Z{P) = Tr = Tr 

\ / Jr'—^OQ 

( 1 ) 

where H = T -I- F is the Hamilton operator. In other 
words, the origin of the method is the notion that the 
finite temperature density matrix , which corre¬ 

sponds to the square of the wavefunction at low and 
to the Maxwell-Boltzmann distribution at high tempera¬ 
ture, can be decomposed into a product of density matri¬ 
ces, each at higher effective temperature Pp = P/P. In 
any case, Eq.[T]is a direct consequence of the Trotter the¬ 
orem and implies that in the limit P —> oo, sampling Zp 
classically is equivalent to the exact canonical quantum 
partition function Z 941. 
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Inserting P — 1 complete sets of position eigenstates 
and introducing momenta using the standard Gaussian 
integral, as well as the symmetric Trotter splitting to 
decompose P, 

Zp = A/-y r y p , (2) 


which can be readily sampled by MD. Herein, = 

is a normalization constant, while r and p are 
the positions and momenta of all N particles. The so- 
called bead-Hamiltonian Pp({r}, {p}) that describes the 
interactions between all iV x P beads, reads as 


Pp({r},{p}) = ^ 


/c=l Li^l 


I rn,ujp / (fc) _ j.(fc+i)' 

(3) 




where P is the number of imaginary-time slices, rrii are 
the particle masses and Wp = P//3 = f3~p^ is the angu¬ 
lar frequency of the harmonic spring potential between 
the beads. As a consequence of the trace of Eq. [U 
Pp({r},{p}) is isomorphic to a classical closed ring- 
polymer, thus [^ . 

For the purpose to reduce the computationally dom¬ 
inating effort to evaluate the intermolecular long-range 
electrostatic interactions P times, the ring-polymer con¬ 
traction scheme is employed [o^. To that extent, the 
Hamiltonian is split into its inter- and intramolecular 
contributions, where the former is limited to a single 
Ewald sum at the centroid of the closed ring-polymer: 


r. 


1 

P 


p 



k=l 


(4) 


At variance to the original PIMD method, the partially 
adiabatic centroid MD (ACMD) approach permits to ap¬ 
proximately compute even dynamical properties within 
the path integral scheme [9g. To that extent, the ef¬ 
fective masses of the ring-polymer beads are chosen so 
as to recover the correct dynamics of the ring-polymer 
centroids. Specifically, the elements of the Parrinello- 
Rahman mass-matrix are selected so that the vibrational 
modes of all ring-polymer beads, except for the centroid, 
are shifted to a frequency of 


n = 


pp/p-1 

f3h 


(5) 


which allows for integration time-steps close to the ionic 
resonance limit [^ . 


Hydrogen Bond Definition and Kinetics 


ria have been proposed 36, Following Skin¬ 

ner and coworkers, we use the potential of mean force 
(PMF) as calculated from the radial-angle joint distribu¬ 
tion function to defi ne a H-bond without any empirical 
parameters [4^, Il0l| . Specifically, the equipotential re¬ 
gion of the 2D PMF surface, which passes through the 
saddle point and encircles the minimum, corresponds to 
the H-bonded state, i.e. the H-bond indicator variable 
B{RoOt0.oho) = where aoHO the angle between 
Roo = |ro,a — ro,d| and the covalent OH bond vector 
(rn — ro)- For all other combinations of Roo and a, 
B{Roo,o;oho) = 0 . 

2 \ For the purpose to study the H-bond kinetics of the 
ipdividual water layers, we define the layer-specific H- 
/bond autocorrelation function 


Cll{t) = {Li{to)Li{to -I- r))/(L*), (6) 

where Li(t) is equal to 1 if at time t a particular molecule 
resides in layer i and 0 otherwise. As a consequence, 
Cll{t) is the correlation that a H-bond, which is existing 
at time to in layer i, is also existing at time to + t within 
the same layer. Thus, Cll{j) eventually relaxes to zero. 


Instantaneous Water/Vapor Interface 

In order to study the structure and in particular the 
dynamics of water as a function of distance from the 
surface requires a reliable definition of the water/vapor 
interface itself. However, using a surface definition 
based on the time-averaged density profile, the rele¬ 
vant spatial fluctuations in space and time of the in¬ 
terface location are neglected. Therefore, the recently 
proposed method to locate the instantaneous interface of 
Willard and Chandler is employed here [t^. Instead of 
a time-averaged density-field, a coarse-grained but time- 
dependent density-field pcg{v,t) is constructed in terms 
of Gaussian functions located at the center of mass of the 
water molecules Vi[t)-. 

Pcg{r,t) = ) 2e A ^ ) , (7) 

where ^ = 2.4 A is a system-dependent coarse-graining 
length. Instantaneous interfaces can now be defined as 
2-dimensional manifolds s(t) = r, for which pcg{s,t) = c. 
In analogy to the Gibbs dividing surface, we have set 
the critical parameter c to equal half of the bulk water 
density. Moreover, it is possible to deduce the proxim¬ 
ity Ui of the i’th water molecule from the instantaneous 
interface for every time-step as 

ai = {[s(t) - r,{t)] ■ n(t)}|s(t)=s.(t)> (8) 


In order to determine whether two water molecules are 
H-bonded or not, different energetic and geometric crite- 


where Si{t) denotes the point on s(t) that is closest 
to while n is the surface normal vector at that 
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point. Averaging over all instantaneous interfaces, the 
ensemble-averaged interface and the corresponding mean 
proximity 


ai = {[(s) -Ti] ■ (n)}|(s>=<s>, 
is obtained. The mean density profile 


p{d) 


1 


N \ 

'^S{ai -d)) , 


(9) 


( 10 ) 


where L is the length of the simulation cell, can be either 
defined in terms of or di to quantify the probability of 
finding a water molecule at a certain distance d from the 
mean or instantaneous water/vapor interface. 

To characterize the orientational dependence of the wa¬ 
ter molecules on the distance to interface, we employed 
the joint conditional distribution: 

P{u,u'\d) = - d)S{cos{e[^'‘) - u) 

X d{cos{eP) - u')'^ , (11) 

where is the angle between n(t) on Si{t) and one of 

the two OH bond vectors r ■ of the i’th 

f 2) 

molecule, while ’ is defined analogously for the other 
OH bond vectors of the same water molecule. High corre¬ 
lations between the surface proximity and the orientation 
of the water molecules appears as peaks in P{u,u'\d). 


COMPUTATIONAL DETAILS 


The AIMD simulation was performed in the canonical 
(NVT) ensemble at 300 K using the second-generation 
Car-Parrinello MD method of Kiihne et al. asimple- 
mented in the CP2K/Quickstep code 68. 102 1. The 

model of the water/vapor interface consisted of a bulk 
water part with 384 light water molecules in a periodic 


s 3 

orthorhombic box of dimension 15.64x15.64x46.92 A . 
To that, an additional vacuum portion of 8 A on both 
sides along the non-periodic z-direction was added, while 
the corresponding Poisson pr oblem was tackled by an ef¬ 
ficient Wavelet-based solver [l03l| . The settings of the 
density functional theory (DPT) calculations to compute 
the interatomic forces were identical with those of our 
previous AIMD studies of water |^, Il04l | . The whole 

system was equilibrated for 1.25 ns by classical MD us¬ 
ing the empirical SPC/Fw interactions potential [s^, 
followed by a DFT-based re-equilibration consisting of 
50 ps, before statistics was eventually accumulated for 
additional 250 ps. The latter was decomposed into 12 
statistically independent 20 ps long trajectories, which 
were separated by 1 ps each. 

The classical MD and quantum mechanical PIMD sim¬ 
ulations were conducted using a flexible water model, 



FIG. 1. Mean density profile of the water slab as obtained 
by the binned time-averaged particle density, as well as from 
the coarse-grained representation pcg(v,t). The dashed line 
denotes the Gibbs dividing surface, which is shown for com¬ 
parison. 


which was parameterized by force-matching to accurate 
DFT-based reference calculations usi^ the TPSS- D3 ex¬ 
change and correlation functional [s^, llOSl Il06l| . The 
corresponding computations were all started from a well 
equilibrated periodic cubic simulation box with II^ water 
molecules each and a density of 0.997 g/cm^. The water 
surface was then created by putting two additional empty 
simulation cells of equal size along the z-direction. The 
resulting system was then re-equilibrated for I ns within 
the NVT ensemble at T=298 K. Thereafter, 200 MD and 
130 ACMD statistical independent 20 ps long trajecto¬ 
ries were computed in the microcanonical (NVE) ensem¬ 
ble. Throughout, periodic boundary conditions were em¬ 
ployed. The short-range interactions were truncated at 
9 A, while the Ewald summation technique was used to 
compute the long-range electrostatic interactions. In case 
of the ACMD simulations, the ring-polymer contraction 
scheme was employed with a cutoff value of 5 A to reduce 
the electrostatic potential energy and force calculations 
to a single Ewald sum, while all other interactions were 
computed using P = 32 beads. The evolution of the 
ring-polymer in time was performed using a discretized 
time-step of 0.1 fs, whereas an integration time-step of 
0.5 fs was employed for the MD and AIMD simulations, 
respectively. 


RESULTS AND DISCUSSION 

Density Profile of the Water/Vapor interface 

In Fig. [T] the time-averaged density profile of the water 
slab as a function of the z-coordinate is shown together 
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Proximity [A] 


FIG. 2. Proximity-based density profiles in terms of the in¬ 
stantaneous ai (top panel) and mean proximity at (bottom 
panel). The instantaneous water layers (L0-L2) are indicated 
by vertical dashed lines. 


with its coarse-grained counterpart pcg{r,t). Based on 
the instantaneous interface definition, the proximity- 
based instantaneous and mean density profiles, are de¬ 
picted in Fig. [2] It is apparent that the impact of NQE 
on the proximity-based density profiles are minor. Fur¬ 
thermore, using the force-matched TPSS-D3 water po¬ 
tential, they are more structured than those of the DFT- 
based AIMD simulation, which exhibits a lower first peak 
and a slightly larger penetration into the vapor phase. 
This may be attributed to the fact that compared to 
bulk water, the average 0-0 distance is increased be¬ 
cause of re laxat ion effects at the water/vapor interface 
'Ml45l.[6i[l07|. Interestingly, most non-polarizable in¬ 


teraction potentials predict the opposite, i.e. a shorten¬ 
ing of the average 0-0 distance. As already described 
by Willard and Chandler, the density profile with re¬ 
spect to the instantaneous interface exhibits pronounced 
peaks within the interfacial region, separated by distinct 
minima [t^. These peaks indicate that close to the in¬ 
stantaneous water/vapor interface, the water molecules 
are arranged in a layered fashion. Therefore, we will re¬ 
fer to the three distinct intervals of 3 A each, as shown 
in the top panel of Fig. [2l as instantaneous water lay¬ 
ers (L0-L2). However, it is important to emphasize that 
LO cannot be understood as a genuine water layer, but 
rather as a sparse population of water molecules with 
a higher proximity to the vapor than to the first water 
layer, which is LI. In fact, despite their equal volumes, 
LO consists only of about 10 — 20% of the water molecules 
in LI. With this novel definition of the instantaneous wa¬ 



FIG. 3. Joint probability distribution as obtained by AIMD 
(left column), ACMD (middle column) and classical MD 
(right column) simulations. The various rows are associated 
with the considered proximity values of d=0.2 A to represent 
LO, 1.2 for Llll and 1.8 A for LI, as well as 4.4 A for L2. 


ter layers in the vicinity of the water/vapor interface it 
is now possible to compute the structure and dynamics 
of interfacial water for each of the various layers individ¬ 
ually, as will be shown in the following. 


Orientation of water molecules close to the surface 

To study the orientation of the water molecules close to 
instantaneous water/vapor interface we employ the joint 
conditional distribution P{u,u'\d). The corresponding 
plots as a function of cos(0) with respect to n(t), are 
shown in Fig. [31 The various rows in Fig. [3] corresponds 
to four distinct proximity values between d=0 A and 5 A 
in order to represent the instantaneous water layers LO- 
L2. All layers beyond L2 do not obey any structural order 
and therefore correspond to bulk liquid water. As before, 
the influence of NQE is minor. Moreover, apart from the 
much enhanced statistics of our MD and PIMD simula¬ 
tions using the force-matched TPSS-D3 water model, no 
appreciable difference to the DFT-based AIMD could be 
detected. 
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In LO (d=0.2 A), cos(0) is either 1 or close to -0.4, 
which corresponds to an angle of 0° and 109°, respec¬ 
tively. From this follows that in the top-most water layer, 
the OH bonds are sticking out of water into the vapor 
phase, but are also pointing into the subjacent layer. In 
this LI (d=1.8 A), the maxima of the joint probability 
distribution are antipodal to LO. Specifically, cos(0) is 
either -1 or approximately -1-0.3, which equates to an an¬ 
gle of 180° and 70°, respectively. As a consequence, in 
LI the OH bonds are on the one hand pointing straight 
into the bulk phase, but on the other hand also towards 
LO, where the water molecules possess lone pair orbitals. 
While this scenario has already been predicted earlier 
41 . 4^ . we now also find an additional region between 
LO and LI where cos(0) « 0, i.e. 0 ~ 90°. We will refer 
to this intermediate layer, which was first predicted by 
Morita and coworkers , where the water molecules are 
preferably parallel oriented with respect to the instanta¬ 
neous water surface as LI 11 (d=1.2 A). Due to the fact 
that both OH bonds are perpendicular with respect to 
the instantaneous surface normal, both of their lone pair 
orbitals are pointing up- and downwards towards LO and 
LI, where the associated water molecules due have their 
OH bonds. Moreover, in LO and LI an additional popu¬ 
lation at cos(0) close to ±0.3 is noticeable that belongs 
to water molecules, where both OH bonds obeys an angle 
of around 70° and 110°, respectively. Put in other words, 
beside water molecules that are parallel with respect to 
the instantaneous water/vapor surface, we also find con¬ 
figurations that are tilted by around ±20°. The latter im¬ 
mediately suggest the possibility of in-plane H-bonding 
between the water molecules of LI 11 and the tilted con¬ 
figurations. The orientation of the water molecules in L2 
is increasingly similar to bulk water and on the onset to 
become fully disordered. Nevertheless, a small orienta¬ 
tional correlation similar to LO is still detectable. 


Hydrogen Bond Network Strnctnre 

In order to investigate the structure of the H-bond net¬ 
work between the water molecules and to eventually de¬ 
vise a model of the water/vapor interface, we compute 
all of our ensemble-averaged H-bond network descriptors 
individually for all of the various instantaneous water 
layers. Specifically, in Table U the average number of 
H-bonds, H-donors and H-acceptors per water molecule 
are given together with the relative fraction of inter- and 
intralayer H-donors, as well as the relative occurrence 
of ND, SD and DD H-bonds. As can be seen, the im¬ 
pact of NQE is again relatively small. Moreover, in the 
bulk water region, the MD simulations employing the 
TPSS-D3 interaction potential and DFT are essentially 
identical. However, in the vicinity of the water/vapor in¬ 
terface, the results between the present water model and 
explicit AIMD simulations differ quantitatively, which is 
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FIG. 4. Layer-specific H-bond autocorrelation function as 
obtained by conventional MD, quantum-mechanical ACMD 
and AIMD simulations. 


a manifestation of the enhanced transferability of elec¬ 
tronic structure based AIMD calculations. 

More importantly, we find that the water molecules 
in LO form on average two H-bonds only. The fact that 
roughly 75% of these H-bonds are SD configurations is 
consistent with our observation that in LO one of the 
OH bonds is protruding out of the water into the vapor 
phase. As already mentioned, this implies that the other 
OH bond is pointing towards LI, which manifests itself 
that more than 80% of the H-bond donors are so-called 
interlayer H-bonds that are connecting LO with LI. 

In LI, the water molecules are forming on average more 
than three H-bonds. Of these, nearly 75% are DD, while 
still approximately 25% are SD configurations. More¬ 
over, at variance to L2 and in particular to LO, about 
70% of the H-bonds are formed in-layer within LI. The 
large fraction of intralayer H-bonds is mainly due to the 
water molecules of LI 11 that are oriented parallel to the 
water surface, while the purpose of the remaining inter¬ 
layer H-bonds is to connect the water molecules of LI 
with those of LO and in particular with L2. This implies 
that the water molecules in LI, which are possessing only 
slightly less H-bonds than in the bulk, are able to com¬ 
pensate the lack of potential H-bond partners owing to 
the presence of the vapor phase by forming additional 
intralayer H-bonds within Lfl/ 

At variance, the H-bond network of the water 
molecules in L2 is rather similar to bulk water with on av¬ 
erage nearly four H-bonds per molecule. Of these, around 
85% are DD configurations. The number of interlayer 
and intralayer H-bonds in L2 is nearly level, with a small 
preference for the latter. 

The H-bond dynamics as determined by Cll{t) is 
shown in Fig. 01 After just 0.5 ps, nearly half of the 
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H-bonds per molecule 
H-donors per molecule 
H-acceptors per molecule 
Interlayer H-bonds in % 
Intralayer H-bonds in % 
Non-donor H-bonds in % 
Single-donor H-bonds in % 
Double-donor H-bonds in % 


(MD/ACMD/AIMD) 

LI (MD/ACMD/AIMD) 

L2 (MD/CMD/AIMD) 

Bulk (MD/ACMD/AII 

2.08 

/ 

2.08 

/ 

1.86 

3.40 

/ 

3.39 

/ 

3.18 

3.74 

/ 

3.73 

/ 

3.68 

3.71 

/ 3.71 

/ 

3.70 

0.96 

/ 

0.96 

/ 

0.84 

1.70 

/ 

1.70 

/ 

1.60 

1.84 

/ 

1.84 

/ 

1.82 

1.84 

/ 1.84 

/ 

1.84 

1.12 

/ 

1.12 

/ 

1.02 

1.70 

/ 

1.69 

/ 

1.58 

1.90 

/ 

1.89 

/ 

1.86 

1.87 

/ 1.87 

/ 

1.86 

0.87 

/ 

0.87 

/ 

0.81 

0.30 

/ 

0.30 

/ 

0.34 

0.46 

/ 

0.46 

/ 

0.44 



/ 


0.13 

/ 

0.13 

/ 

0.19 

0.70 

/ 

0.70 

/ 

0.66 

0.55 

/ 

0.55 

/ 

0.57 


-/- 

/- 


0.13 

/ 

0.13 

/ 

0.20 

0.02 

/ 

0.02 

/ 

0.02 

0.01 

/ 

0.01 

/ 

0.01 

0.01 

/ 0.01 

/ 

0.01 

0.78 

/ 

0.78 

/ 

0.75 

0.26 

/ 

0.26 

/ 

0.35 

0.14 

/ 

0.14 

/ 

0.16 

0.14 

/ 0.14 

/ 

0.15 

0.09 

/ 

0.09 

/ 

0.04 

0.72 

/ 

0.72 

/ 

0.63 

0.85 

/ 

0.85 

/ 

0.83 

0.85 

/ 0.85 

/ 

0.85 


TABLE I. Ensemble-averaged H-bond distribution of the individual interfacial water layers by conventional MD, quantum- 
mechanical ACMD and AIMD simulations. 



FIG. 5. Schematic of the water/vapor interface from the var¬ 
ious MD simulations. The topmost water molecule represent 
LO, while the lowest constitute LI. The interjacent molecules 
correspond to LI'L 


water molecules left LO, while more than 80% of the 
molecules remained in LI and L2, respectively. Hence, 
LO is particularly mobile, whereas the H-bonds between 
the water molecules in LI and especially LI 11 are par¬ 
ticularly strong. Interestingly, employing the TPSS-D3 
water model, the decay of Cll{t) for LO and L2 is more 
pronounced than in our AIMD simulations, while for LI 
this behavior is reversed. This implies that in the AIMD 
simulations, the enhanced H-bond strength of the water 
molecules in LI 11 is less pronounced. 


CONCLUSION 


Altogether, based on our simulations the following 
model of the water/vapor interface, which is depicted in 
Fig. [SJ is eventually emerging: The topmost layer of the 
water surface is dominated by SD water configurations, 
where the dangling OH bonds are preferably sticking out 
of the water into the vapor phase, while at the same time 
serving as an H-bond acceptor and donor for LI. We find 
that this second water layer can be further divided into 
the previously proposed LI 4ll 4^, where the orientation 
of the water molecules is inverted with respect to LO, and 
a novel interjacent layer denoted as LI^L Specifically, in 
the former, one of the OH bonds is preferentially pointing 
towards the bulk-like L2, while at the same time accept¬ 
ing and donating H-bonds from LO and LI 11. In this latter 
water layer, the molecules are oriented parallel with re¬ 
spect to the water/vapor interface and are able to form 
H-bonds with LO and LI, as well as particularly strong in¬ 
tralayer H-bonds within LlH. The water molecules in L2 
are already structurally rather disordered and in general 
resembles bulk water with a relatively weak orientational 
correlation that is similar to LO. All of this implies that 
only the topmost ~ 5 A obeys structural oder. 
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